function [gbest,gbestval,fitcount,fitness] = FCO_func(fhd,Dimension,Particle_Number,Max_Gen,lowerbound,upperbound,pos,bool_Plot,varargin)

gbest = 0;
gbestval = 0;
fitcount = 0;
rand('state',sum(100*clock));
me = Max_Gen;
ps = Particle_Number;
D = Dimension;
g = 0.13;
G0 = 0.28;
epsilon = 1e-3;
cc=[2 2];   %acceleration constants
iwt=0.9-(1:me).*(0.5./me);

num_FClevel = 5;
prob_forage = 0.9;
prob_invdir = 0.5;
dir_forage = 1;
size_each_level = repmat(floor(ps/num_FClevel),1,num_FClevel);
size_each_level(1) = ps - (size_each_level(1)*(num_FClevel-1));
lb_level(1,1) = 1;
ub_level(1,1) = size_each_level(1);
for l = 2 : num_FClevel
    lb_level(l,1) = sum(size_each_level(1,1:l-1))+1;
    ub_level(l,1) = sum(size_each_level(1,1:l));
end

[pub, plb, vub, vlb] = SetBound(D,ps,lowerbound,upperbound);
% pos = VRmin+(VRmax-VRmin).*rand(ps,D);

fitcount = ps;
vel = zeros(ps,D); % initialize the velocity of the particles
acc = zeros(ps,D); % initialize the velocity of the particles

fitness = feval(fhd,pos',varargin{:});
[fitbest,gbestindex] = min(fitness);
fitworst = max(fitness);
gbest = pos(gbestindex,:);
pbest=pos;
pbestval=fitness;
gbestval = fitbest;
gbestrep=repmat(gbest,ps,1);
if(bool_Plot)
    figure(8)
    haxes = plot( 0 , 0 );
    XArray = [1];
    YArray = [gbestval];
    title('FCO');
end
figure(8)
haxes = plot( 0 , 0 );
XArray = [];
YArray = [];
title('FCO');
% level_std = zeros(me,num_FClevel,D);
% return
for iter = 2:me    
    if(fitbest == fitworst)
        mass = repmat(1/ps,ps,1);
    else
        mass = (fitness'-fitworst)/(fitbest-fitworst);
    end
    mass = mass / sum(mass);
    sort_mass = [mass (1:ps)'];
    sort_mass = sortrows(sort_mass,1);
    std_mass = std(mass);
    if(iter>30 & mod(iter,10)==0)
        XArray = [ XArray iter]; 
        YArray = [ YArray std_mass];
        set( haxes , 'XData' , XArray , 'YData' , YArray );
        drawnow
    end
    
    a=2-iter*((2)/me);
    % Decomposer ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
    for l = 1 : num_FClevel
        N = size_each_level(l);
        mean_pos = mean(pos(sort_mass(lb_level(l):ub_level(l),2),:));
        std_pos = std(pos(sort_mass(lb_level(l):ub_level(l),2),:));
%         level_std(iter, l,:) = std_pos;
        level_mass = sort_mass(lb_level(l):ub_level(l),1);
        for n = 1 : N
            index_n = n + lb_level(l) - 1;
            index_n = sort_mass(index_n,2);
            index_perturb = fix(rand * ps) + 1;
            if(l == 1)
                vel(index_n,:) = rand(1,D).*(pub(index_n,:)-plb(index_n,:))+plb(index_n,:);
            elseif(l == 2)
                vel(index_n,:) = mean_pos + (3+rand(1,D)).*std_pos;
            elseif(l == 3)
                vel(index_n,:) = ForageQGSA(g, mean_pos, pos(index_n,:));
            elseif(l == 4)
                index_target = RouletteWheel(level_mass,1);
                index_target = index_target + lb_level(l) - 1;
                index_target = sort_mass(index_target,2);
                vel(index_n,:) = pos(index_n,:) + (pos(index_target,:)-pos(index_n,:)).*(2*rand(1,D)) + ...
                                                  (rand < 0.01)*(pos(index_perturb,:)-pos(index_n,:)).*(8*rand(1,D));
            elseif(l == 5)
                S = rand;
                index_target = (S>0.1)*ub_level(l) + (S<=0.1)*ub_level(l-1);
                index_target = sort_mass(index_target,2);
                vel(index_n,:) = pos(index_n,:) + (pos(index_target,:)-pos(index_n,:)).*(2*rand(1,D)) + ...
                                                  (rand < 0.01)*(pos(index_perturb,:)-pos(index_n,:)).*(4*rand(1,D));
            elseif(l == 6)
                vel(index_n,:) = ForageQGSA(g, pos(sort_mass(ps, 2),:), pos(index_n,:));
                vel(index_n,:) = pos(index_n,:) + (pos(sort_mass(ps, 2),:)-pos(index_n,:)).*(2*rand(1,D));
                vel(index_n,:) = pos(index_n,:) + (pos(sort_mass(ps, 2),:)-pos(index_n,:)).*(2*rand(1,D)) + ...
                                                  ((l~=num_FClevel) | (n~=N))*(rand < 0.01)*(pos(index_perturb,:)-pos(index_n,:)).*(a*rand(1,D));
                vel(index_n,:) = pos(index_n,:) + (pos(sort_mass(ps, 2),:)-pos(index_n,:)).*(2*rand(1,D)) + ...
                                                  (rand < 0.01)*(pos(index_perturb,:)-pos(index_n,:)).*(a*rand(1,D));

%             end
        end
    end    
 
%     vel = ((vel>=vlb)&(vel<=vub)).*vel+(vel<vlb).*vlb+(vel>vub).*vub;
    pos = vel;
    pos = ((pos>=plb)&(pos<=pub)).*pos+(pos<plb).*plb+(pos>pub).*pub;
    
    fitness = feval(fhd,pos',varargin{:});
    sort_fitness = [fitness' (1:ps)'];
    sort_fitness = sortrows(sort_fitness, 1);
    fitbest = sort_fitness(1,1);
    fitworst = sort_fitness(ps,1);
    gbestindex = sort_fitness(1,2);
    gbest = pos(gbestindex,:);
    gbestval = fitbest;
    tmp=(pbestval<fitness);
    temp=repmat(tmp',1,D);
    pbest=temp.*pbest+(1-temp).*pos;
    pbestval=tmp.*pbestval+(1-tmp).*e;%update the pbest
    gbestrep=repmat(gbest,ps,1);
    
    if(bool_Plot)
        XArray = [ XArray iter]; 
        YArray = [ YArray gbestval];
        set( haxes , 'XData' , XArray , 'YData' , YArray );
        drawnow
    end
end

% for d = 1 : D
%     figure(d)
%     for l = 1 : num_FClevel    
%         subplot(2,3,l)
%         plot(level_std(:,l,d))
%         title(['Level ' num2str(l)])
%     end
% end

end
